probit.m

This script calculates population values of Bresnahan and Reiss's estimates of

based on their ordered probit model of entry.

Contents

Undertake the prerequisite calculations, if necessary.

The required inputs are CN, the matrix of possible (C,N) pairs and v, the vector with each pair's probability in the ergodic distribution.

if ~exist('probitARG','var')
    bellman
    thresholds
    ergodic
    CN=[kron(ones(maxN-minN+1,1),omega') kron((minN:1:maxN)',ones(nCstates,1))];
else
    CN=probitARG.cn;
    v=probitARG.v;
end

CData = exp(CN(:,1));
NData = CN(:,2);
Current plot held
Current plot held
Current plot held
Current plot held

Create an in-line function for the standard normal c.d.f.

Phi = @(x) 0.5*(x>=0).*(erf(abs(x)/sqrt(2))+1)+0.5*(x<0).*erfc(abs(x)/sqrt(2));

Eliminate all sates with near-zero probability.

This avoids numerical errors arising from multiplying very, very small probabilities by ``infinite'' log likelihood evaluations.

eps=1e-10;
CData=CData(v>eps); NData=NData(v>eps); v=v(v>eps); v=v/sum(v);

Set up the binary likelihood function and its derivatives.

We first estimate the entry thresholds with standard binary probits. These estimates then serve as starting values for the ordered probit estimation. So that it can be used for all of the thresholds, the rank under consideration is read as one of the likelihood function's inputs.

lf = @(x) v'*((NData>=x(3)).*log(Phi(log(CData/x(1))/x(2))) + (NData<x(3)).*log(Phi(log(x(1)./CData)/x(2))));

% Calculate the minimum and maximum values of N consistent with the ergodic distribution.
Nlow=min(NData);
Nhigh=max(NData);

Maximize the binary likelihood functions.

This uses the optimization toolbox for its minimization.

bProbitC=zeros(Nhigh-Nlow-1,1);
bProbitSigma=bProbitC;
for R=Nlow+1:Nhigh;
    f = @(x) -lf([x R]);
    x=fminsearch(f,[0.5*(exp(overlineC(R))+exp(underlineC(R))) 0.2]);
    bProbitC(R)=x(1);
    bProbitSigma(R)=x(2);
end

Set up the full multinomial likelihood function.

As written, x must be a row vector with [sigma epsilon Cbar(Nlow) Cbar(Nlow+1) ... Cbar(Nhigh) infinity].

logprobs = @(x) log(Phi(log(CData./x(NData-Nlow+2)')/x(1))...
    -Phi(log(CData./x(NData-Nlow+3)')/x(1)));

%Likelihood function. First argument is standard deviation, and remaining arguments are thresholds. (``zero'' and
%``infinity'' coming first and last).
lf = @(x) v'*logprobs([ x(1) 1e-10 x(2:end) 1e10]);

Maximize the ordered probit likelihood function.

f= @(x) -lf(x);
xm=fminsearch(f,[mean(bProbitSigma) bProbitC']);
mProbitC=xm(2:end);
mProbitSigma=xm(1);
 
Exiting: Maximum number of function evaluations has been exceeded
         - increase MaxFunEvals option.
         Current function value: Inf 

Report results if this is a test run.

if ~exist('probitARG','var')
   disp(mProbitC);
   disp(mProbitSigma);
   disp(f(xm));
end
         0    0.7164    1.1383    2.0959

    0.3631

   Inf